#include "defines.h"


int test_z_blocking(DATA_TYPE * B, int NX, int NY, int NZ, int T, int Bx, int By, int Bz, int tb) {

	if(Bx<(2*XSLOPE+1) || By<(2*YSLOPE+1) || Bx>NX || By>NY || tb>min(((Bx-1)/2)/XSLOPE,((By-1)/2)/YSLOPE)){
		return 0;
	}

	DATA_TYPE (*A)[NX+2*XSTART][NY+2*YSTART][NZ+2*ZSTART] = (DATA_TYPE (*)[NX+2*XSTART][NY+2*YSTART][NZ+2*ZSTART]) B;

	int bx = Bx-2*(tb*XSLOPE);
	int by = By-2*(tb*YSLOPE);

	int ix = Bx+bx;
	int iy = By+by;

	int xnb0 = ceild(NX,ix);
	int ynb0 = ceild(NY,iy);
	int xnb11 = ceild(NX-ix/2+1,ix) + 1;
	int ynb12 = ceild(NY-iy/2+1,iy) + 1;
	int xnb12 = xnb0;
	int ynb11 = ynb0;
	int xnb2 = max(xnb11,xnb0);
	int ynb2 = max(ynb12,ynb0);

	int nb1[2] = {xnb12 * ynb12, xnb11 * ynb11};
	int nb02[2] = {xnb2 * ynb2, xnb0 * ynb0};
	int xnb1[2] = {xnb12, xnb11};
	int xnb02[2] = {xnb2, xnb0};

	int xleft02[2] = {XSTART-bx, XSTART+(Bx-bx)/2};
	int ybottom02[2] = {YSTART-by, YSTART+(By-by)/2};
	int xleft11[2] = {XSTART+(Bx-bx)/2, XSTART - bx};
	int ybottom11[2] = {YSTART-(By+by)/2, YSTART};
	int xleft12[2] = {XSTART-(Bx+bx)/2, XSTART};
	int ybottom12[2] = {YSTART+(By-by)/2, YSTART-by};

	int level = 1;
	int tt,n;
	int t, x, y, z;
	int zz;
	int ymin, ymax;
	int xmin, xmax;
	int zmin, zmax;



	for(tt =- tb; tt < T; tt += tb){
#pragma omp parallel for  private(xmin,xmax,ymin,ymax,t,x,y,zz, zmin, zmax)
		for(n = 0; n < nb02[level]; n++){

			for(zz = 0; zz - 2 * tb * ZSLOPE < NZ + ZSLOPE; zz += Bz) {

				for(t = max(tt,0); t < min(tt + 2*tb, T); t++) { 

					xmin = max(   XSTART,   xleft02[level] + (n%xnb02[level]) * ix      - tb*XSLOPE + myabs(t+1,tt+tb) * XSLOPE);
					xmax = min(NX+XSTART,   xleft02[level] + (n%xnb02[level]) * ix + bx + tb*XSLOPE - myabs(t+1,tt+tb) * XSLOPE);
					ymin = max(   YSTART, ybottom02[level] + (n/xnb02[level]) * iy      - tb*YSLOPE + myabs(t+1,tt+tb) * YSLOPE);
					ymax = min(NY+YSTART, ybottom02[level] + (n/xnb02[level]) * iy + by + tb*YSLOPE - myabs(t+1,tt+tb) * YSLOPE);
					zmin = max(ZSTART,            zz      - (t-max(tt,0))*ZSLOPE);
					zmax = min(NZ+ZSTART,          zz + Bz - (t-max(tt,0))*ZSLOPE);
					for(x = xmin; x < xmax; x++) {
						for(y = ymin; y < ymax; y++) {
							#pragma ivdep
							#pragma vector always
							for (z = zmin; z < zmax; z++) {
								kernel(A);
							}
						}
					}
				}
			}
		}

#pragma omp parallel for private(xmin,xmax,ymin,ymax,t,x,y,zz, zmin, zmax)
		for(n = 0; n < nb1[0] + nb1[1]; n++) {
			for(zz = 0; zz - 2 * tb * ZSLOPE < NZ + ZSLOPE; zz += Bz) {
				for(t = tt + tb; t < min(tt + 2*tb, T); t++) {
					if(n < nb1[level]) {
						xmin = max(     XSTART,   xleft11[level] + (n%xnb1[level]) * ix       - (t+1-tt-tb) * XSLOPE);
						xmax = min(NX + XSTART,   xleft11[level] + (n%xnb1[level]) * ix  + bx + (t+1-tt-tb) * XSLOPE);
						ymin = max(     YSTART, ybottom11[level] + (n/xnb1[level]) * iy       + (t+1-tt-tb) * YSLOPE);
						ymax = min(NY + YSTART, ybottom11[level] + (n/xnb1[level]) * iy  + By - (t+1-tt-tb) * YSLOPE);
					}
					else {

						xmin = max(     XSTART,   xleft12[level] + ((n-nb1[level])%xnb1[1-level]) * ix      + (t+1-tt-tb) * XSLOPE);
						xmax = min(NX + XSTART,   xleft12[level] + ((n-nb1[level])%xnb1[1-level]) * ix + Bx - (t+1-tt-tb) * XSLOPE);
						ymin = max(     YSTART, ybottom12[level] + ((n-nb1[level])/xnb1[1-level]) * iy      - (t+1-tt-tb) * YSLOPE);
						ymax = min(NY + YSTART, ybottom12[level] + ((n-nb1[level])/xnb1[1-level]) * iy + by + (t+1-tt-tb) * YSLOPE);

					}
					zmin = max(ZSTART,            zz      - (t-(tt + tb))*ZSLOPE);
					zmax = min(NZ+ZSTART,          zz + Bz - (t-(tt + tb))*ZSLOPE);
					for(x = xmin; x < xmax; x++) {
						for(y = ymin; y < ymax; y++) {
							#pragma ivdep
							#pragma vector always
							for (z = zmin; z < zmax; z++) {
								kernel(A);
							}
						}
					}
				}
			}
		}
		level = 1 - level;
	}
}
